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SECTION  1 


INTRODUCTION 

In  chin  report  w«  prnnnnc  «  general  aodnl  for  the  notions  of  a 
tvo-dlaana Iona 1  airfoil  in  unsteady  aerodynaaics.  The  equations  of 
notion  are  a  sec  of  coupled  ordinary  and  functional  differential 
equations.  Making  use  of  general  "approximating  Vagner  functions" 
we  develop  a  state  space  aodel  that  can  be  used  for  identification 
and  optlaal  control  of  auch  ayateaa. 

In  Section  2  we  formulate  the  basic  (hereditary)  equations  and 
In  Section  3  we  develop  the  state  apace  model.  Although  the  resulting 
model  is  infinite  dimensional,  there  are  certain  special  features  of 
the  state  space  model  that  can  be  exploited  In  order  to  construct 
approximation  schemes.  In  Section  4  we  discuss  these  features  and 
present  an  approximation  technique  that  can  be  used  for  parameter 
estimation  and  optimal  control  of  such  systems.  Section  5  contains 
numerical  examples  that  illustrates  the  basic  ideas.  Section  6  presents 
some  preliminary  numerical  results  from  an  analysis  of  data  derived 
from  wind-tunnel  tests  on  an  oscillating  airfoil.  Throughout  this 
report  we  shall  use  the  symbol 

dB  f  (a)  g  (s) 

to  denote  the  Rlemann-Stleltjes  Integral  of  g  with  respect  to  f  over 
the  Interval  [a,b]  (see  [6]  for  a  dlseuasion  of  Rlemann-Stleltjes 
integration) . 
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SECTION  2 


THE  HEREDITARY  MODEL 

Consider  the  two-dlasasional  airfoil  shown  la  figure  1.1*  The 
basic  equations  for  the  pitching  and  plunging  notions  of  the  air¬ 
foil  aay  be  found  in  a  nuaber  of  papers  (see  for  example  [10,17])* 
We  shall  follow  the  development  in  [17]  and  use  fairly  standard 
nomenclature  and  conventions.  We  shall  not  write  down  the  exact 
constants  which  may  also  be  found  in  the  various  references  cited 
previously. 


Let  h(t)  denote  the  "plunge"  end  a(t)  the  "pitch"  of  the  air¬ 
foil  at  tine  t.  The  equations  of  notion  can  be  written  in  the 
form 

M#  z(t)  +  K#  z(t)  -  -F(t),  (1) 

where  s(t)  is  the  coluan  vector  z(t)  ■  col[h(t)  o(t)]  end  F  con¬ 
tains  the  aerodynaaic  loads  on  the  airfoil.  In  particular. 


where  L  end  M  are  the  aerodynaaic  loads  corresponding  to  total 
a 

wing  lift  per  unit  depth  end  total  aoaent  about  the  1/4  chord  per 
unit  depth,  respectively.  For  the  airfoil  considered  here,  it 
follows  that  (note  that  in  this  case  d  •  ■  1/2) 

Ma(t)  -  vpb3  {[j  i(t>  +  [0  O]  i(t)>,  (3) 

end 

L(t)  -  rpb2  {[1  £]  ii(t)  +  [0  0]  s(t>)  (4) 

+  (2*pUb)  D(t), 

where  D(t)  is  the  "Duhenel  integral."  In  particular,  D(t)  is 
given  by 

D(t)  -  /*  4(g(t  -  t))  Q(x)dT,  (5) 


Mmvmia i  u  .i.i  ui  i  u  ■  . hi  mjmiij.i  i.imii  i  iu.uui.mi  wmMp 

.  .  . . ^  _  . 


where  4(Ut/b)  Is  the  Wagntti  function,  and 

Q(t)  -  /  {[1  bl  i(t)  +  [0  0]  *(T)}dT.  (6) 

To  obtain  a  state  space  model  that  is  suitable  for 
identification  and  optimal  control,  one  must  provide  a  "useful" 
representation  for  D(t).  One  approach  to  this  problem  is  to  ap¬ 
proximate  D(t)  by  approximating  the  Vagner  function  4.  For 
example,  the  Jones  two-term  exponential  approximation  of  the 
Wagner  function  (see  [13,14})  results  in  a  finite  dimensional 
model  where  D(t)  may  be  viewed  as  the  output  of  a  second  order 
linear  control  system.  Recently,  it  has  bean  shown  that  D(t) 
may  be  realized  4s  the  output  of  an  infinite  dimensional  control 
system  (see  [1,  2,  10]).  In  particular,  D(t)  can  be  explicitly 
represented  as  the  output  of  a  system  governed  by  a  functional 
differential  equation  (l.e.,  a  hereditary  system,  see  [9]).  If 
this  explicit  representation  is  exploited,  then  the  resulting 
model  becomes  an  infinite  dimensional  control  system.  Again,  it 
will  be  necessary  to  Introduce  some  type  of  approximation  in  order 
to  obtain  a  workable  finite  dimensional  model. 

Although  we  shall  make  use  of  the  fact  that  D(t)  is  the  out¬ 
put  of  an  hereditary  system,  the  approach  presented  below  is  sim¬ 
ilar  in  Api/uct  to  the  Jones  two-term  exponential  approximation 
of  the  Wegner  function  which  we  now  briefly  describe. 

Jones  [13]  approximated  the  Wagner  function  4  by  the  two- 
term  exponential  functlor 
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(7) 


n4.  ~0.  Ut/b  -0,  Ut/b 

1(f)  -  1  -  ^  e  1  -y  2 

where  "  0.165,  02  ■  0.335,  0^  ■  0.0455  end  02  "  0.3.  If  4  is 
substituted  into  (5) ,  then  D(t)  is  approximated  by 


D(t)  -  QCt)  -  B1(t)  -  a2  B2<t), 

where  B^,  B2  satisfy  the  ordinary  differential  equations 

r 

®l(t)  "  "  $1  b  Bl(t)  +  ^(t) 

< 

B2(t)  -  -  02  2.  B2<t)  +  Q(t)  . 


(8) 


(9) 


This  approach  leads  to  a  sixth  order  ordinary  differential  equation 
model  for  the  motion  of  the  airfoil  (see  [17]  for  complete  details). 

In  order  to  understand  the  hereditary  model  to  be  discussed 
below,  it  is  worthwhile  to  look  at  the  same  approximation  scheme 
(i.e.,  the  approximating  model  resulting  from  equations  (7) -(9)) 
from  a  different  (system  theory)  point  of  view.  Vrlte  the  Wagner 
function  as 

*(£t)  -l-W(t),  (10) 


and  observe  that  the  Duhamel  integral  (5)  is  equal  to 

D(t)  -  Q(t)  -  y(t)  (11) 


where 


y(t)  -  /  W(t  -  t)  Q(T)dT  . 
0 


(12) 


The  function  W(t)  Is  celled  e  u&igkting  pcubtMn  end  the  Jones  two- 
term  exponential  approximation  (7)  Is  equivalent  to  approxi¬ 
mating  W(t)  by 


U/b 


W(t)  -  [a^  o2l  e 


h 0 1 

i°  -»2j  ri 


H 


(13) 


Observe  that  (13)  can  be  written  as 


W(t)  -  C  X(t)  B  , 


where  X(t)  Is  the  fundamental  matrix  for  the  ordinary  differential 
equation  (of  dimension  2) 


(14) 


C  -  [a^  a^]  and  B  -  colli  ll.  In  particular. 


D(t)  -  Q(t)  -  y(t) 


(15) 


where  y(t)  is  the  output  of  the  second  order  control  system 

i(t)  -  b  [  o1  ♦[i]«(t) 


(16) 


with  output 
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(17) 


y(t)  -  [<*!  a2I  *(*> 

Consequently,  one  aay  view  the  Jones  approxlaetion  as  an  approxi¬ 
mating  of 

y(t)  -  /*  W(t  -  t)  Q(t)<It 
0 

by  the  output  of  the  (second  order)  finite  dimensional  system 
(16)-(17)  with  input  $(t). 

In  the  model  detailed  below,  we  propose  to  epproxlmate  y(t) 
by  the  output  of  a  simple  hereditary  system  (l.e.,  a  retarded 
functional  differential  equation) .  There  are  two  basic  reasons 
for  considering  this  approach.  As  indicated  above,  one  can  show 
that  the  model  should  include  hereditary  terms  (see  [2],  [9]). 
Moreover,  since  ordinary  differential  systems  may  be  viewed  as 
"special"  functional  differential  equations,  the  inclusion  of 
delayed  terms  is  a  very  natural  extension  of  the  "exponential 
type"  approximations. 

He  turn  now  to  the  general  problem  of  approximating  y(t) 

(and  hence  D(t))  by  the  output  of  a  simple  hereditary  system.  For 
each  parameter  p  €  R11,  let  G(t,p)  denote  a  n  *  n  matrix  valued 
function  such  that  for  p  6  Ru  the  mapping  t  •*  G(t,p)  is  of 
bounded  variation  on  compact  subintervals  of  (-•,0J.  Let  B  and 
C  be  n-dimenslonal  colon  and  row  vectors,  respectively.  Let 
r  >  0  and  consider  the  n-dlmensional  retarded  functional  differ¬ 
ential  equation 
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(18) 


0 

*d(t)  ■  I  dgG(«,p)xd(t  +  s)  +  BQ(t). 

The  output  of  system  (18)  Is  defined  by 

yd(t)  -  C  xd(t).  Q9) 

The  system  defined  by  (18)-(19)  is  parameterized  by  y  -  (p,r, 

T 

B,C  ).  We  define  a  to  be  the  set  of  all  parameters  y  “  (p,r,B, 

CT)  €  x  K  x  S.n  x  Kn  such  that  r  >  0  and  the  pair  B,  C 
satisfies 

n 

CB  -  l  c.b.  «  1/2  .  (20) 

1-1  1  1 

The  set  0  is  called  the  set  of  admiaible.  paJvamJtXM  for  (18)- 
(19). 

Given  (p,r),  let  X(t;p,r)  denote  the  n  x  n  matrix  satisfying 

(  •  o 

X(t)  -  /  d  G(s,p)X(t  +  s)ds,  t  >  0, 

-r  * 

X(0)  -  In,  (21) 

X(s)  -  0n,  s  <  0, 

\ 

vhere  I  and  0  denote  the  n  x  n  identity  and  zero  matrices, 
n  n 

respectively.  The  matrix  X(t;p,r)  Is  called  the  {,wndamejvCcU. 
matxix  for  (18).  (see  [7,  12  J).  Given  y  -  (p,r,B,CT)  6  0, 
define  the  approximate  UagneA  function  «(Ut/b:y)  by 
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(22) 


-  1  -  W(t;Y), 

where  W(t;y)  la  tba  weighting  pattern  for  the  hereditary  aye tea 

(18)  -(19)  given  by 

W(t;y)  ■  CX(t;p,r)B.  (23) 

Obaarva  that  the  condition  that  CB  ■  1/2  iapllea  that 
•(DO/bjy)  ■  1/2.  Horeover,  if  y(t;Y)  denotes  tha  output  of  (18)- 

(19) ,  then  it  follows  that  (aee  Hale  [12]) 

y(t;Y)  -  /*  «(t  -  t;y)  Q(t)dt.  (24) 

0 

Let  D(tJY)  and  L(t;Y)  he  defined  by  (aee  equatlona  (4)-(S)) 

D(t;y)  •  /  *(r(t  -  t))  Q(t)dt  ■  Q(t)  -  y(t;y)  (25) 

0  0 

and 

L(C;y)  -  »Pb2  {[1  |]  a(t)  +  [0  0]  i(t)>  (26) 

+  (2vpUb)  D(t;Y>. 

respectively.  If  L(f,y)  la  substituted  for  L(t)  in  (l)-(2), 
then  we  obtain  the  nodal 

\  ii(t)  +  a(t)  -  -F(t;Y)  (27) 


where 
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SECTION  3 


THE  STATE  SPACE  MODEL 
T 

Given  the  paraaetar  y  -  (p,r,B,C  )  €  Q  va  shall  construct  a 
dynaalcal  systaa  for  tha  airfoil  described  in  Figure  1,1,  Bacall 
that  s(t)  -  col[h(t)  a(t)]  and  *d(t)  is  the  n-dlnsnslonal  vec¬ 
tor  function  that  satisfies  tha  retarded  functional  differential ; 
equation  (18).  Let  x(t)  denote  tha  (4  +  n) —dimensional  vector 

x(t)  -  col[  6(t)  o(t)  h(t)  a(t)  x.  (t) ] .  (29) 

U 

The  basic  nodal  is  defined  by  ((27) -(28)) 

Mg  a(t)  +  Kb  i(t)  -  -F(t;y)  (30) 

where 

Db  “(t>  <31> 

Observe  that 


D(t;y) 


T 

(Q(t)  -  y(t;y)} 

LoJ 


and  under  tha  assuaptlon  that  s(0)  -  s(0)  -  0,  It  follows  that 


[ol?<tiT>*[o  )]*»*['•  I]*<t)'[l],<t!T)-  <5l) 

Using  ths  fsct  thst  y(t;r)  ■  C  *d(t),  it  follows  tram  aquations 
(31)  and  (32)  that  (30)  Is  aquivalant  to 

f°  ”]i 

Lo  Ub J 


M  *(t)  +  K  *(t)  -  -irpb2  I  2  |  s(t)  -  tpb2 
•  •  1  b_  3b_ 

2  8 


-  2irpUb 


[o  o]  i(l>  -  J”“  [1  I] 

[  0  ]  X^(t)  • 

Dafina  ths  2x2  aatrlcas  Ell*  *11  *12  by 


s(t)  (33) 


s(t) 


+  2wpUb 


“11 


r,  b  i 

+  *pb2 

1  2 

■ 

b  3o 

L  2  8  J 

(34) 


*11  “  ~ 


_ 

2 

0  U 

wpo 

0  Ub 

■  «* 

+  2wpUb 


1  b 
0  0 


(35) 


and 


*12’- 


+  2*pUb 


n 


(36) 


raspactlvaly.  If  danotaa  tha  2  *  n  aatrlx 


12 


-  2irp0b 


[el-*~[*  C2"'C-1 

[o  J  L°  0  ...  0  J 


(37) 


then  the  syetea  (33)  can  be  written  ee 


*11  “  *ii  *(*>  +  ^  *<e>  +  ^3  *d(0< 


(38) 


Moreover,  Xj(t)  eetiefiee 


id(t)  -  /_r  dfl6(e,p)zd(t  +  a)  +  BQ  (t) 


(39) 


-  /_r  d#G(e,p)*d(t  +  s)  4-  B{[1  b]ic(t)  +  (0  0]i(t)}. 


Define  the  n  *  2  eetrlcee  E^  end  A^  by 


e31  -  -B(l  b)  -  - 


r  bx  bb3  1 
*>2  bb2  ■  • 


b  bb  -« 
n  n 


(40) 


end 


A^  -  B{0  0)  - 


ro  Obj 

0  Ob- 


0  Ob  J 

D 


(41) 


reapectlvely.  With  end  A^  eo  defined,  equation  (39)  beeoaea 


E31  *(t)  +  *d(t)  -  dgG(e,p)*d(t  +  s)  ♦  An  i(t),  (42) 
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(43) 


la  tana  of  tha  "state  ▼actor"  x(t)  -  col[z(t)  z(t)  *d(t)J, 
aquations  (38)  and  (42)  combine  to  fon  tha  systaa 

0 

Ex  (t)  -  A*  (t)  +  /  d  H  (a,p)  x  (t  +  a), 

-r  ■ 

where  E,  A  and  H(s,p)  are  defined  by 


(44) 


(45) 


(46) 

Obaarve  that  la  general  (43)  la  a  linear  retarded  functional 
differential  aquation  with  flalta  delay  r  >  0.  lowawar,  the  nodal 
that  raaulta  from  tha  Janes  two-tan  exponential  apprewineClen  la 
but  a  special  ease  of  (43) . 

BNfU  3.1.  Lot  a  ■  2,  p  •  (-•. »  -lj)  and  C  ■  (a^  Sj].  fix 
r  >  1  and  let  >  *  sail  1  D .  Safina  tha  2  »  2  aatrln  wnlnad  funa- 
tlan  «(a,p)  hp 


E11 

0 

0 

E  - 

0 

*2 

0 

e_. 

0 

I 

31 

n 

\l 

*12 

*13 

A  - 

X2 

0 

0 

*31 

0 

0 

and 


B(s,p) 


0  0  0 

0  0  0 

0  0  G(a,p) 
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if  •  -  0, 


if  ■  <  0  . 


If  Xj(t)  is  say  continuous  2-dimensional  vector  function 
on  C-r,  +») ,  then 


0  _ 

t  d  G(s,p)x  (t  +  s)  - 
-r  8  a  b 

end  the  "delay  equation"  (39)  reduces  to  the  ordinary  differential 
equation 


+1  ° 


0  -6, 


*d(t>. 


*d(t>  *  ? 


■®1  0 

0  -8, 


*d(t)  + 


Q(t). 


The  output  in  this  case  is  given  by 


y(t;y)  -  t«x  «2,xd(t)* 


and  the  corresponding  fundamental  matrix  X  (t)  -  X (t;y)  is  given  by 


X(t;y)  - 


-UAb  Bxt 


-U/b  82t 


Consequently ,  the  approximating  Wagner  function  defined  by  (20)-(23) 
ie  given  by 
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_  -D/b  fLt 

!•  1  o 

[al  °2]  [  0  ,’0/b  V 

-U/b  St  -0/b  B,t 

<*!  •  -  «2  *  » 

which  is  the  Jones  two-tens  exponential  approximation .  Therefore, 
the  general  model  (43)  includes  as  special  cases  those  models 
derived  from  two-term  "Jones  type"  approximations. 

EXAMPLE  3.2.  Assume  that  the  airfoil  represents  a  tail  section 
of  a  flexible  aircraft  and  that  this  tail  section  is  in  the  wake 
of  the  forward  wing.  The  aerodynamic  loads  on  the  tail  at  time 
t  are  influenced  by  the  motion  of  the  forward  wing  at  some  pre¬ 
vious  time,  say  (t-r).  Moreover,  because  of  the  elastic  coupling 
between  the  tall  and  the  wing,  the  notion  of  the  forward  wing  is 
influenced  by  the  loads  on  the  tail.  Consequently,  the  loads  on 
the  tall  at  time  t  Involve  the  motions  of  the  tail  at  a  previous 
time  (t-r).  This  "time  delay"  feature  of  a  wing-tall  configuration 
has  been  observed  by  Reding  and  Ericsson  [  14,15]  in  their  stability 
analysis  of  the  747/Orblter. 

To  see  how  one  can  Incorporate  these  pure  time  delays 
in  the  model  defined  by  (3.15)  we  assume  that  r  >  0  and  n  >.  1  are 
fixed.  Let  B  and  C  be  such  that  CB  *  1/2  and  for  p  €  assume 
that  Aq(p),  A^(p)  are  n  *  n  matrices  parameterized  by  p.  Define 
the  n  x  n  matrix  function  G(s,p)  by 


*(£t)  -  *(£t  ;y)  -  1  - 
-  1  - 


l 

! 
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If  x^(t)  la  a  continuous  n-dlmensional  vaetor  function,  than 
0 

/  dgG(a,p)*d(t  +  a)  -  AQ(p)xd(t)  +  A1(p)xd(t  -  r). 

Consequently,  tha  functional  differential  aquation  (18)  becomes 
the  delay-differential  equation 

Xj(t)  -  AQ(p)xd(t)  +  A1(p)xd(t  -  r)  +  B  Q  (t) . 

Also,  the  basic  model  (43)  is  equivalent 

Ex(t)  -  FQ(p)x(t)  +  Fx(p)  x(t  -  r) 

where  E  is  defined  by  (44)  and  Fg(p)  and  F^(p)  are  given  by 


*11  *12  *13 

F0(p)  -  I2  0  0 

L*31  0  *0<P>. 

and 

’  0  0  0 

Fj(p)  -  0  0  0  , 

0  0  A1(p) 


(48) 


(49) 


respectively 


We  conclude  this  section  with  a  model  that  includes  a  pur*  time 
delay  aa  wall  aa  distributed  delays.  This  aodal  Includes  all  of  the 
previously  discussed  aodals  and  yet  a  maker  of  nice  features  are 
retained  which  nake  the  nodal  a  useful  tool  for  Identification  and 
optimal  control. 

Let  r  >  0  and  n  >  lba  fixed  and  CB  •  1/2  be  as  above.  Assume 

that  6  (s.p)'la  an  n  x  a  matrix  function  continuous  la  a  and  p  and 
c 

define  G(s,p)  by 

G(e,p)  -  Cg(a,p)  +  /_*  Ge(u.p)du, 
where  Gg(s,p)  Is  defined  by  (47),  l.e. 

• 

Ap(p)  +  Aj(p),  a  -  0, 

Ga(s,p)  -  Ajfo),  -r  <  a  <  0,  (50) 

0,  s  <_  -r. 

If  C(s,p)  la  defined  as  above,  and  x^(t)  Is  a  continuous  n-dinsn- 
alonal  vector  function,  then 

I  °  dgC(a,p)xd(t  +  s)  -  AQ(p)zd(t) 

0 

+  A1(p)xd(t  -  r)  +  /  Gc(o,p)xd(t  +  a)ds. 

Consequently,  equation  (43)  nay  be  written  as 

Sx  (t)  ■  Fq( p)x(t)  +  rx(p)x(t  -  r)  ♦  /  K(s,p)x(t  +  s)ds,  (51) 


where  £  is  defined  by  (44),  FQ(p),  F^(p)  ere  defined  by  (48)-(49) 
and  K(stp)  is  given  by 

ro  o  o  1 

K(e,p)  -  0  0  °  (52) 

Lo  0  Gc(s,p)J 

Although  the  model  described  by  (51)  has  the  nice  feature 
that  it  in  some  sense  generalizes  the  classical  approach  due  to 
Jones  [13],  there  are  also  theoretical  justifications  for  consid¬ 
ering  (51)  as  the  basic  model.  In  particular,  the  "delayed 
state"  Xj(t)  is  closely  related  to  the  circulation  function  I*(t) 

(see  [4,9])  and  it  may  be  shown  using  theoretical  aerodynamics 
that  r(t)  does  in  fact  satisfy  a  linear  functional  differ¬ 
ential  equation.  A  complete  discussion  of  these  theoretical 
aspects  will  appear  in  a  forthcoming  paper  (see  [9]), 

The  structure  of  the  system  (51)  is  very  special.  There 
are  four  "non-delayed  states"  (h(t)  ,  a(t),  h(t)  ,  a(t))  and  n 

"delayed  states"  x£j(t) .  More  precisely,  the  matrices  F^(p)  and 
K(s,p)  consist  of  zeros  except  in  the  lower  right-hand  n  x  n 
blocks.  This  particular  feature  can  be  exploited  in  the  develop¬ 
ment  of  numerical  algorithms,  resulting  in  enormous  savings  in  storage 
and  CPU  time.  In  the  next  section  we  present  two  numerical  schemes 
for  Identification  and  control  of  systems  governed  by  equations  of 
the  form  (51) .  These  schemes  are  modified  versions  of  the  so 
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called  AVE  and  SPLIHE  schemes  (aae  [3,4,5])  chat  have  been 
developed  especially  for  general  hereditary  system  of  the  form 
(51). 


r 
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SECTION  4 

THE  REDUCED  STATE  MODEL  AND  NUMERICAL  APPROXIMATIONS 

In  this  section  we  consider  the  hereditary  control  system 
governed  by  the  linear  functional  differential  of  order  (4  -I-  n) 

E(p)x(t)  -  FQ(p)x(t)  +  Fx(p)x(t  -  r)  (33). 

0 

/  K(s,p)x(t  +  s)ds  +  Gu(t) , 


with  initial  data 


x(0)  -  n,  x(s)  -  $(s),  s  <  0,  (54) 


and  output 


y(t)  *  Hx(t), 


(55) 


where  n  is  a  (4  +  n) -dimensional  vector  and  $(t)  is  a  (4  +  n)- 
dimensional  vector  valued  function.  Moreover,  we  assume  that 
the  system  matrices  E(p),  F^(p)  and  K(s,p)  have  the  special 
forms. 


S(p) 


En(p)  0  0 

0  i2  0 
e3i(P)  0  in 


F1(p) 


0 

0 

0 


0  0 

0  0 

o  A1(p) 


(56) 


(57) 
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and 


K(e,p)  ■ 


0 

0 

0 


0 

0 


0 

0 


0  Gc(a,p) 


(58) 


We  assume  that  for  each  p  €  »U  the  n  *  n  matrix  valued  function 
s  •+•  Gc(s,p)  is  integrable  on  [-r(0]  and  the  2x2  matrix  Eu(p) 
is  non-singular. 

Since  Ej^ (p)  is  non-singular,  it  follows  that  E(p)  is  non¬ 
singular  and 


[E^f1  - 


'  [En(p)]"1  0  0 

0  I2  0 

-e31(p)[e11(p)]‘1  o  in 


(59) 


The  basic  idea  la  to  approximate  the  hereditary  system 
(53) -(55)  by  a  (possibly  large)  system  of  ordinery  differ¬ 
ential  equations.  This  approach  to  Identification  and  optimal 
control  of  hereditary  systems  haa  been  used  by  a  number  of  in¬ 
vestigators  (see  [3]  for  a  survey  of  this  technique). 

There  are  two  particular  schemas  that  have  been  used  for  such 
hereditary  systems;  the  Al/E  scheme  (a  "finite  difference"  type 
scheme  to  be  described  below)  and  the  SPLINE  schema  (a  "finite 
element”  type  scheme,  see  [4,5]). 
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Although  thM*  two  srhamas  h«v«  boon  applied  to  o  nuabor  of 
hereditary  systeas,  there  la  a  problem  with  the  procedure.  In 
particular,  the  also  of  the  approximating  ays tea  of  ordinary  dif¬ 
ferential  equatlona  can  becoaa  very  large  thus  reducing 
effectiveness  of  the  algorlthaa.  However,  the  special  fora  of 
the  hereditary  systea  (53)-(55)  can  be  exploited  to  greatly  im¬ 
prove  the  usefulness  of  these  general  algorithms.  To 
describe  the  necessary  aodif lcatlons ,  we  shall  concentrate  on  the 
scheae.  Similar  modifications  are  valid  for  the  SPLINE  scheme 
and  a  complete  theoretical  and  numerical  analysis  of  both  proce¬ 
dures  (Including  convergence  raaulta)  will  appear  In  a  future 
paper. 

Let  0  <  ra  <  rM  <  4«  and  suppose  that  F  cl11  x  la 

a  compact  convex  set.  He  assume  that  the  Initial  data  n  6 
and  [-r^.O]  ♦  are  fixed  and  consider  the  following 

typical  parameter  estimation  problem  (see  [8]). 

problem  (ID).  Given  the  input  function  u  and  obeeAvations 
yi  at  tune*  (0  <  ^  <  F2  <  .  .  .  <  t^  <  T) ,  find  the  &y&tejn 
panametesu  y*  -  (p*,r*)  6  r  which  minimize  the  tfiAOK 

l  II  7<ttiy)  -  y.  ||2,  (60) 

i-1  1  x 

MfceAe  y(t;y)  ii  the  output  at  time  t  to  the  kcuditaxy  eyetem 

(53)-(55). 


23 


Given  y  -  (p,r)  €  r  and  a  poaitive  integer  M,  partition  tha 
intarval  [-r,0]  into  H  aubiataxvala  [t*(r) ,  t*_j(r)],  where 
tj (r)  ■  "Jt/H,  J  B  0|  1|  2| •  • •  i  !•  Tor  J  ■  0t  1>  2» .  .  »  •  M 
define  +*  by 

-  n.  (61) 

and 


6 


M 

j 


u  t"  .(r) 

“  t  J-* 

r  N 

tj(r) 


4(a)ds. 


(62) 


Note  that  4"  la  an  (4  -f  n) -dineaaional  vector  and  that 


H _ .  .  ,N  .M  .H. 

,  col  [*  ,  4  *  • •  •  *  4  I 
J  ’  0  i  H 


(63) 


la  an  (»  +  1)(*  +  n)-dinenaional  vector. 

Define  the  l(H  +  1)(4  +  n)J  *  [(N  +  1)(4  +  n)]  aatrlcea 
EN(p)  and  AN(p)  by 


Ap) 


E(p)  0 
0  0 


(64) 


and 
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typ>  kJ(p>  kJ(p)  •••  ^(p)  P1(p)+rJJ(p) 


—i  z!ij  o 

r  r  u 


A"(p)  -I  0  Ji 


where 


\(P)  m  J  *  K(s,p)ds. 
J  J*/  \ 

'j<r) 


L*t  H  and  G  be  the  (H  +  1)(4  +  n) -dimensional  row  and  coluan 


vectors  defined  by 


h"  -  [H,  0,  0,  ....  0] 


G  -  cel  [0,  0,  0,  ,  0], 


respectively. 


With  «*(p),  AH<p),  xj,  H*  and  G*  defined  by  (61)-(68) 
shove,  the  **  apptotimatoig  contort  intern  is  defined  by 

lH(p)iM(t)  -  AH(p)xH(t)  +  G*u(t) 
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vlcb  Initial  condition 


n 


i 

\ 

t 


and  output 


*  (0> 


(70) 


yK(t)  -  H“x"(t). 


(71) 


Corraapondlng  to  system  (69) -(71)  we  consider  the  following 
approximating  parameter  estimating  problem. 

PROBLEM  (IDN) .  Seven  the  input  4/unction  u  and  obieevatiom 
7t  at  time*  t±  (0<  t,  <  t2  <  ...  <  tm  <  T) ,  4ind  the.  iyitem 
pafumeteei  y*  -  (pN,  rH)  €  r  which  minimize  the  vuuok 


Ay)  -  |  I  II  J*6±iy)  -  y.  ||2.  (72) 

2  i-1  1  1 

wheAe  yM(t;y)  is  the  output  at  time  t  to  the  W**1  approximating 
intern  (69) -(7i). 

Under  very  reasonable  assumptions  it  is  possible  to  show 

*N 

that  u  N  •*•  »,  the  optimal  solutions  y  to  PROBLEM  (IDS)  con¬ 
verge  to  y*.  the  optimal  solution  of  PROBLEM  (ID)  (see  [3.8]). 
Therefore,  the  idea  is  to  select  a  value  of  N.  construct  the 
approximating  system  (69)-(71)  and  use  a  standard  algorithm 
to  solve  PROBLEM  (IDN)  for  yS.  If  N  is  "large  enough",  then  y* 
is  a  good  approximation  to  y*.  This  particular  numerical  scheme 
Is  called  the  AVE  scheme. 
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Consider  the  sis*  of  chs  approximating  system  (69) -(71)  de- 
finsd  by  tbs  Al/E  scheme  above.  The  dimension  of  this  systsa  is 


DDi(M,n)  -  (N  +  1) (4  +  n) ,  which  for  large  values  of  H  (or  n)  grows 
very  rapidly.  This  is  illuatrated  in  Table  4.1  for  various  values 
of  N  and  for  n  ■  1,  2,  3. 


TABLE  4.1.  The  diaenaion  of  the  AVE  approximating  syataas 
for  values  of  N  and  n. 


N 

DIM(N.2) 

DIMQI.3) 

2 

15 

18 

21 

4 

25 

30 

35 

8 

45 

56 

63 

16 

85 

102 

119 

32 

165 

198 

231 

64 

325 

390 

455 

For  N  -  16,  32,  64  ...  ,  the  slse  of  the  approximating  system 
(69) -(71)  becomes  an  important  (and  limiting)  consideration  in 
numerical  coaputatlons .  Therefore,  we  seek  soae  method  to  reduce 
this  computational  burdon.  To  modify  the  Ai/E  schema  it  is  sufficient 
to  make  a  simple  but  important  observation  concerning  the  special 
structure  of  the  hereditary  equation  (53).  In  particular,  we  note 
that  at  moat  n  of  the  (4  +  a)  states  in  (53)  are  delayed.  Consequently 
it  is  not  necessary  to  approximate  all  of  the  "history"  x(t  -  r), 
but  only  those  coord Inst as  which  are  delayed. 
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-HUssm: 


Lac  Bj  danot*  the  marl  mu.  ^mber  of  noa-saro  oalane  la 

Aj(»)  and  the  n«bar  of  non-ter,  c  ms  In  6e(a»p).  The  BMfcar  n4 

aaac  fall  kaewaan  0  and  n.  If  *  J,  than  A^(p)  ■  0  and  c  p)  I 
0  ao  chat  (S3)  la  an  ordinary  diffarantlal  aquation  and  no  approx 
laaflaa  ia  naadad.  Conaaquantly,  wa  aaaoaa  that  1  <  <  n  «n<i 

partition  tha  aatrlcaa  FQ(p) ,  F^(p)  and  K(e,p)  aa  follows : 

yp>  •  t»01(p)  »  VM<p».  (73) 

»j(p)  -  [0  i  Fu(p)l.  (74) 

K(a,p)  •  (0  !  Kj(a,p)),  (73) 

'*i»  7n(p),  YU<P>  hj(a.p)  ar*  1(4  ♦  n)  ■  n4J  natriaaa. 

■an  n  lava  aaanaad  at  Chant  laaa  af  panarality  that  tha  fine 

(4  ♦  a)-u.  aelaaaa  an  tha  aero  ealana.  Lat  [Oil.)  daaaca 
*  *d 

tha  a4  »  (4  ♦  a)  antrix  renal  acini  af  the  aatrix  with  aana  la 

tha  first  (4  a  n)-n4  oalanaa  and  dafino  p”,  j  •  0,  1,  2,  .  .  .  ,  ■ 

hy 

♦J  -  n  (74) 


[0  )  I 

“d 


a7) 
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•hero  ♦*  1*  defined  by  (62).  Hoc*  that 

wjj  -  col  [#*,  (78) 

la  m  [«d*l  +  (4  +  a)  ]-d1— nolonal  vector. 

Define  tho  [n^  +  (4  +  «)]  *  [n^-H  +  (4  +  #)]  aatrlcos 

«"(»)  and  aJ<»)  by 

‘  E(p)  0  0  .  .  . 

0  I  0  .  .  . 

m  Ui 

eJ(p)  -  .  in  .  .  . 

•  •  d 

o  o  t 

0 

0  0  .  .  . 

and 


0 


0 


(81) 


where  I  represents  Che  n .  x  a.  Identity  end 


Kjj  *  /  J  K2(*,p)de. 
t*  (r) 


The  fact  that  only  n^  of  the  states  In  the  hereditary  equation 
(53)  are  delayed  laplies  that  the  state  space  model  defined  by  (53)- 
(55)  Is  not  a  "minimal  realization"  for  the  system.  Consequently, 
the  corresponding  Hth  approximating  system  defined  by  (69) -(72) 

Is  not  a  minimal  (*.a.  not  controllable  and  observable)  and  can  be 
reduced  In  dimension.  In  particular.  If  we  define  the  REDUCED  W**1 
appJWXAjnating  iyitvn  by  the  ordinary  differential  equation 


^(p)vN(t)  -^(p)wN(t)  +  G*u(t)  ,  (82) 


with  Initial  data 


N,_.  ** 

w  (0)  -  Wq 


(83) 


and  output 


yj(t)  -  HjAt)  ,  (84) 

where  hJ  -  [H,  0 . 0]  and  G*  -  col  [G,  0,  .  .  .  ,  0],  than  one 

can  show  that  (69) -(71)  and  (82) -(84)  have  the  same  Input- 
output  operator  (l.e.  y*(t)  -  yN(t)). 
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The  dimension  of  the  reduced  systea  defined  by  (82)-(84)  depends 
on  N,  n  end  n^.  Moreover,  the  ayetea  ia  largest  when  n^  -  n.  For 

applicatlone  to  the  airfoil  in  Figure  1.1,  the  value  of  n^  ■  1  ia 

/ 

reasonable  and  in  this  case  the  dimension  of  the  reduced  systea  (82)- 
(84)  la  given  by  DXM(N,n)  ■  [N  +  (4  +  n) ] .  Table  4.2  illustrates* 
the  sizes  of  the  REDUCED  appKoxJjncuUng  iy-itvni  for  various  values  of 
N  and  n,  in  the  case  that  n^  ■  1.  If  one  compares  Table  4.2  with 
Table  4.1,  then  it  becoaea  clear  that  for  large  values  of  N  and  n 
the  systea  (82) -(84)  is  much  smaller  (and  hence  much  more  manageable 
for  numerical  calculations)  than  the  non-reduced  AVE  system  (69) -(72) 
For  example,  if  H  »  64  and  n  ■  3  then  the  dimension  of  (69)-(72)  is 
455  while  the  dimension  of  the  "equivalent”  system  (82) -(84)  is  71! 

TABLE  4.2.  The  dlasnslon  of  the  REDUCED  AVE 
approximating  systems  for  values 
of  N  and  n. 


nd“  1 

N 

DIMR(N.l) 

DIMR(N,2) 

DIMR(N.3) 

2 

7 

8 

9 

4 

9 

10 

11 

8 

13 

14 

15 

16 

21 

22 

23 

32 

37 

38 

39 

64 

69 

70 

71 

3: 


Corresponding  to  problems  (ID)  and  (IDN)  ve  have; 

PROBLEM  (IDNr)  .  Given  thz  input  function  u  and  obAZAvatixms 
y.  at  timzA  t.  (0  <  7.  <!,<...  <t  -  T) ,  £ind  tint  &y*tzm  pa*a*> 

1  1  —  1  Z  B 

am  aw  am 

metZAA  yr  ”  (pR,  rR)  €  r  tohich  minimize,  thz  zMjok 

jJ(y)  -  j  II  Yr^jy)  -  y±  ||2,  (85) 

tohzfie.  yR(t;Y)  i&  thz  output  at  timz  t  to  thz  fuzdaczd  W^1  approx¬ 
imating  Aystzm  (82)-(84). 

Although  we  restricted  ourselves  to  the  AVE  approximating 
scheme,  the  same  basic  ideas  can  be  applied  to  the  SPLINE  scheme. 
However,  there  are  some  differences  in  the  two  algorithms.  For 
example,  it  is  not  true  that  the  REDUCED  SPLINE  SYSTEM  has  the 
same  input-output  operator  as  the  SPLINE  SYSTEM. 


SUCTION  5 


NUMERICAL  EXAMPLES 


In  this  section  ve  present  some  numerical  examples  to  Illus¬ 
trate  the  method  discussed  above.  All  the  numerical  results  were 
run  under  the  Conversational  Monitor  System  on  an  IBM  370/158  at 
the  Virginia  Tech  Computer  Center.  The  test  model  is  described 
by  the  second  order  delay  equation 


.  rx.(t)i  o  lnTx/oi  ro  oi  K(t-r)i  m 

4  1  -  ,  1  ♦  1  + 

I_x-(t)j  l_-u  oj  Lx,(t)J  |_0  -aj  l_x,(t-r)J  L 1 J 


with  initial  data 


^  1  .  PI. 

L  v>  J  L 0  J 


a  <  0 


(87) 


and  output 


*  TjCO  " 

'  1  0  ' 

— 1 

X 

/—s 

rr 

_ 1 

.  ?2(t)  - 

_  0  1  . 

L  x2(t)  J 

(08) 


where  u_^(t)  is  the  unit  step  at  t  ■  0.1.  The  final  time  of 
T  *  2  is  used  in  each  example. 

The  "observations"  y^  were  generated  by  selecting  the  "true" 
set  of  parameters  y*  "  (u* ,  a* ,  r*)  -  (4 ,  10 ,  1) ,  using  the  method 


(86) 
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of  steps  to  solve  (86)-(88)  for  y(t)  and  settin'  y,  ■  y  (t.,). 
Using  this  procedure,  we  generated  101  equally  spaced  data  points 
on  the  Interval  [0,2]. 

For  various  values  of  N  ■  2,  4,  8,  16,  .  .  .  the  two  problems 
PROBLEM  (IDN)  and  PROBLEM  (IDNR)  were  solved.  The  amount  of  CPU 
time  required  to  solve  each  problem  was  recorded  in  order  to  com¬ 
pare  the  two  algorithms.  In  addition  to  the  AVE  scheme  detailed 
above,  we  tested  the  SPLIHE  scheme  so  that  a  comparison  of  these 
two  schemes  could  also  be  made.  The  maximum  likelihood  algorithm 
described  in  [  8 ]  was  used  In  all  of  the  calculations. 
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EXAMPLE  5.1 

In  this  example  we  assumed  that  the  time  delay  r  ■  1  Is  known 
and  estimate  u  and  a.  Since  the  maximum  likelihood  estimator  Is 
an  Iterative  algorithm,  start  up  values  have  to  be  provided .  For 
this  example  the  initial  estimates  were  u  »  /IF  and  a  ■  8,0.  The  nu¬ 
merical  results  for  this  example  are  summarized  In  Tables  (5.1)  - 
(5.2). 

Observe  that  the  AVE  scheme  did  not  converge  for  N  ■  2,  4. 

For  N  -  32  the  dimension  of  PROBLEM  (IDN)  is  66.  Consequently, 
(because  of  storage  limitations)  this  approximating  system  Is 
"two  large"  and  we  could  not  solve  PROBLEM  (IDN)  for  N  -  32.  The 
AVE  scheme  produces  exactly  the  same  parameter  estimates  as  the 
REDUCED  AVE  acheme.  However,  the  REDUCED  AVE  scheme  requires  con¬ 
siderably  less  CPU  time  than  the  AVE  scheme. 

The  SPLINE  algorithm  converged  for  all  values  of  N.  Moreover, 
as  expected  the  SPLINE  scheme  produced  better  parameter  estimates 
than  did  the  AVE  scheme.  Note  that  not  only  does  the  reduced  prob¬ 
lem  require  less  CPU  time,  the  REDUCED  SPLINE  scheme  produces  bet¬ 
ter  parameter  estimates  than  the  full  SPLINE  scheme.  Figures  5.2 
through  5.3  Illustrate  typical  "converged"  data  fits  for  N  ■  8, 

16,  32  using  the  AVE  scheme.  Figure  5.4  shows  the 
N  ■  8  SPLINE  data  fits.  The  figures  are  plots  constructed  by  using 
the  reduced  approximating  systems.  However,  it  Is  Impossible  to 
distinguish  these  plots  from  the  plots  one  obtains  using  the  noe- 
reduced  approximating  systems  snd  therefore  those  plots  are  not 
Included. 
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TABLE  5.1  A  comparison  of  parameter  estimates  obtained  using 
the  AVE  and  REDUCED  AVE  approximations. 


PROBLEM  (IDN)  - 

AVE 

PROBLEM  (IDNR)  - 

AVE 

N 

-N 
a> _ 

-N 

a 

CPU 

N 

*N 

*N 

CPU 

2 

DID 

NOT 

CONVERGE 

2 

DID 

NOT 

CONVERGE 

4 

DID 

NOT 

CONVERGE 

4 

DID 

NOT 

CONVERGE 

8 

3.2476 

15.5113 

47  sec 

8 

3.2476 

15.5113 

20  sec 

16 

32 

3.6431 

12.2373 

171  sec 

16 

32 

3.6431 

3.8233 

12.2373 

11.1275 

47  sec 

189  sec 

True 

4.0000 

10.0000 

True 

4.0000 

10.0000 

TABLE  5.2 

A  comparislon  of  the  parameter  estimates  obtained  using 
the  SPLINE  and  REDUCED  SPLINE  annroxlmatloris. 

PROBLEM  (IDN)  - 

SPLINE 

PROBLEM  (IDNr)  - 

SPLINE 

N 

-N 

03 _ 

-N 

a 

CPU 

N 

*N 

-N 

fR 

CPU 

2 

3.7738 

7.8675 

18 

2 

3.6466 

9.1217 

8 

4 

3.9674 

9.9133 

17  sec 

4 

3.9535 

10.3377 

13  sec 

8 

3.9981 

9.9651 

36  sec 

8 

3.9950 

10.0643 

20  sec 

16 

32 

3.9979 

9.9731 

94  sac 

16 

32 

3.9970 

3.9976 

10.0000 

9.9842 

41  sec 

184  sec 

True 

4.0000 

10.0000 

True 

4.0000 

10.0000 
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G1R2N16R 


I  TR=  9 


FIGURE  5.2.  Th«  data  fits  for  Examp la  5.1  using  N  »  16  and  the  AVE 
approximations . 


T  (2) 


01R2N16R  I  T R  =  9 


FIGURE  5.2.b.Tha  data  fits  for  Exaopia  5.1  uaing  N  *  16  and  tha  AVE 
approxlnatlona . 


0 1 R2N32 A 


I  TR=  9 


T  (2) 


G 1 R2N32R  I  TR=  9 


FIGURE  5.3.b.Thc  data  fita  for  Example  5.1  ualng  N  -  32  and  Che  AVE 
approximations . 


T  (2) 


0 1 R2N8SP 


I  TR=  8 


FIGURE  5 . 4 J> .  The  data  fits  for  Example  5.1  using  N  »  8  and  the  SPLINE 
approximations . 


EXAMPLE  5.2 


V*  now  consider  the  problea  of  estlaating  all  throe  paraaetera 
(w(  a,  r)  in  the  aodel  (86)-(88).  Start  up  values  were  taken  to 
be  u  -  Sl5,  a  -  8.0  and  r  -  0.8.  The  addition  of  the  tlae  delay  as 
an  unknown  paraaeter  Increases  the  complexity  of  the  problem  and  the 
AVE  scheme  never  conversed.  However,  the  SPLINE  scheae  converged 
for  all  adalsalbla  values  of  R.  Tables  5.3  and  5.4  contain  a 
auaaary  of  the  nuaerlcal  results  for  this  example.  As  In  the  pre¬ 
vious  example,  the  reduced  problea  required  less  CPU  tlae  to  con¬ 
verge  and  generally  produced  better  paraaeter  estlaates  than  the 
full  SPLINE  scheae.  Figures  5. 5  through  5.  7  Illustrate  the 
converged  data  fits  for  N  •  2,  8,  32. 
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TABLE  5.3.  Example  5.2  parameter  estimates  using  SPLINE. 


SPLINE 

SCHEME  -  PROBLEM  (IDN) 

N 

-N 

«N 

*N 

r 

CPU 

fa) 

a 

2 

3.4340 

8.1849 

.8725 

31  sec 

4 

3. 8313 

9.0540 

.9476 

28  sec 

8 

3.9949 

9.9439 

.9989 

60  sec 

16 

3.9986 

9.9775 

1.0002 

241  sec 

True 

4.0000 

10.0000 

1.0000 

TABLE  5.4 .  Example  5.2  parameter  using  REDUCED  SPLINE. 
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FIGURE  5.6.b.The  data  fits  for  Example  5.2  using  N  ■  8  and  Che  SPLINE 
approximations . 
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SECTION  6 


PRELIMINARY  ANALYSIS  OF  WIND-TUNNEL  DATA 

Several  numerical  experiments  war*  conducted  using  actual  wind- 
tunnel  data  for  a  system  alallar  to  that  shown  in  Figure  1.1.  The 
experiments  wars  performad  at  Nellsen  Engineering  and  Research,  Inc. 
In  Mountain  View,  CA,  under  the  supervision  of  Dr.  S.  C.  McIntosh. 
Dr.  McIntosh  gratioualy  made  soma  of  tha  data  available  to  the 
authors. 

The  mathematical  modal  appropriate  for  the  experimental  system 
differs  slightly  from  that  analysed  above.  Flretly,  tha  elastic 
axis,  at  which  point  the  plunging  notion  la  measured,  is  not  at  the 
quarter-chord  point  (xa  t  1/2).  Thus,  while  the  ness  matrix  Ma  has 
the  fora 


Ma- 

_  S«  i«_ 

the  aerodynamic  added  mass-matrix  la  somewhat  more  general  than  the 
form  used  In  (3.6) 


irpb2 


1  b  •  x 

a 

b  •  xa  b2(l/S  +  x2) 


A  second  generalization  is  the  inclusion  of  viscous  damping  tens 

and  CQ  along  the  diagonal  in  the  submatrix  A^  (see  (35)).  The 

structural  stiffness  matrix  is  diagonal  with  elements  k  and  k  ,  as 

n  a 

in  [17]. 

Since  the  elastic  axis  is  not  necessarily  at  the  quarter-chord, 
the  aerodynamic  moment  must  include  a  term  involving  the  product  of 
the  lift  and  a  moment  arm  b(xQ  -  1/2).  The  aerodynamic  variable  now 
"feeds  into"  the  pitch,  as  well  as  the  plunging  equation  of  motion. 
Thus  the  A^  sub-matrix  is  written  (see  (37)). 


A13  -  2ir  p  U  b 


C. 


The  generalization  to  non-coincident  elastic  axis  and  aerodynamic 
center  also  requires  alteration  of  the  sub-matrix  which  now  has 
the  form  (see  (40)) 


E31  -  -  B[l,  b(l/2  +  xa)]  . 

The  nominal  values  of  the  various  system  parameters  are  given 
in  Table  6.1 


54 


TABLE  6.1.  Nominal  values  of  the  system  parameters, 
parsmeter  value 

1.628* 10"3lb-sec2/in2 
5 . 426 •10_3lb-sec2 
6.741* 10“* lb-sec2 /in 
4.95  In 
.297 

1.47  *  10“7lb-sec2/in* 

6.944  lbs/in2 
43.1  lbs 

1.9* 10-3lb-aec/ in2 
8.7*10~3lb-sec 
1325  In  In/sec. 


Note  that  the  only  nodel  parameters  for  which  values  have  not  been 
specified  are  the  matrices  B  and  C-and  the  matrix-valued  function 
G(s;p).  These  quantities  describe  the  dynamic  model  for  the  unsteady 
aerodynamic  behavior  (see  (21)-(23)>. 
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EXAMPLE  6.1 


The  first  exenple  for  which  we  present  results  employs  e  first- 
order  ordinary  differential  equation  model  for  the  unsteady  aero¬ 
dynamic  behavior.  This  amounts  to  using  a  single  exponential  to 
approximate  the  weighting  pattern  W(t)  in  the  Wagner  function  repre¬ 
sentation  2.10.  Thus,  we  take  n  ■  1,  B  -  1  and  G(s;p)  ■  g^u^s), 
where  is  the  usual  unit  step  at  zero.  In  this  case  the  Riemann- 
Stleltjes  integral  in  (42)  becomes 

r  <i*l81u0(s)]xd(t+s)  -  gj*d(t)  . 

so  that  (42)  is  an  ordinary  differential  equation.  Note  that  from 
(20)  theoretically  one  should  expect  C  ■  1/2.  With  this  information 
the  overall  system  (43)  is  a  fifth-order  ordinary  differential 
equation. 

The  data  used  consists  of  sampled  values  of  the  variables  h(t), 
a(t),  h(t)  and  o(t)  at  157  equally  spaced  points  on  the  Interval  [0, 

.25  sec.].  Since  the  data  seemed  relatively  free  of  noise,  the 
Initial  values  for  the  first  four  components  of  the  state  were  taken 
as  the  recorded  values  at  the  initial  time.  No  initial  value  is 
available  for  state  component  x^ ,  the  aerodynamic  variable,  thus  it 
was  identified.  Note  from  (35)  that  the  second  matrix  in  the  sum 

for  includes  a  multiplicative  factor,  2ir,  which  is  theoretical  i 

1 

lift-curve  slope.  In  this  experiment  the  lift-curve  slope  was 

i 

i 

< 


Identified  ao  chat  the  (1,1)  and  (1,2)  elements  la  the  sub-matrix 

A^  are  allowed  to  wary  accordingly.  The  other  parameters  Identified 
are  the  scalar  coefficient  C,  which  appeara  la  the  aub -matrix  Aj^, 
and  the  aerodynamic  time  constant  g^. 

Insults  of  this  sxparlasat  ara  llluatrated  In  Table  6.2  sad  In 
Figures  6.1  and  6.2.  Zn  these  figures  y(l)»li,  y(2)«a,  y(3)-h  and 
y(4)  -  ou 

TABLE  6.2.  Parameter  estimates  for  tbs  ordinary  dlfferenltal 


aquation  nodal;  Example  6.1. 


Iteration 

Xj(0) 

C% 

C 

*1 

0 

0.0 

6.28 

.500 

-12.2 

1 

25.8 

11.65 

.838 

-28.6 

5 

95.6 

8.61 

.726 

-95.2 

10 

86.3 

8.68 

.652 

-75.6 

20 

78.9 

8.98 

.628 

-66.0 

Shown  la  Figure  6.1  are  the  data  matchea  ualng  the  start-up 
estimates  for  the  parsmatsrs.  As  Is  evident  from  the  y(l)  and  y(2) 
matches,  the  fundamental  frequency  Is  nearly  correct  but  the  actual 
plunga  damping  Is  higher  then  predicted  by  the  modal  while  the  pitch 
dmsplng  Is  lighter.  Tbs  actual  data  exhibits  some  high  frequency 
oscillations  that  are  completely  missing  la  the  modal.  The  "converged" 
results  shown  la  Figure  6.2  Indicate  sons  improvement  In  the  data  fits 
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but  the  quality  of  the  fit  is  still  quite  poor.  It  seeas  likely  that 
soae  of  the  "fixed  parameters"  (e.g.  the  plunge  damping  parameter,  C^) 
may  be  in  error.  Alternatively,  it  may  be  that  the  assumed  model 
structure  (aerodynamics  via  a  first-order  ordinary  differential  equation) 
is  not  appropriate.  In  any  case  a  more  systematic  investigation  is 
needed.  As  in  the  previous  figures  we  use  "  +  +  +  "  to  denote  the 
wind-tunnel  data  values  and  a  solid  line  to  indicate  the  values  of  the 
output  from  the  mathematical  models. 
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T  (3) 


SfiM3 


c.  The  startup  data  fits  for  the  ordinary  ditfert 


0.00 

FIGURE  6.1. 


T  (2) 


2.50 


1.50 


0.50 


-0.50 


-1.50 


-2.50 


0.00 

FIGURE 


r  (3) 


SM  1  3  I  T R  =  10 


FIGURE  6.2.c.  The  converged  data  fits  for  the  ordinary  differential  equation 
model;  Example  6.1. 
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EXAMPLE  6.2 


A  second  sxpsrlasnt  wee  conducted  using  s  sleple  daisy  model  for 
the  unsteady  aerodynamic  behavior.  The  situation  Is  very  similar  to 
the  ordinary  differential  equation  case  described  above,  except  that 
the  matrix  valued  function  G(s;p)  now  has  the  form  G(s;p)  ■  g^Ug(s) 

+  g2u_r(s)  so  that  the  Riemann-Stlalt jes  Integral  In  equation  (42): 
now  becomes 


/°r  ds  G(s;p)xd(t+s)  -  gjXjCt)  +  gjX^Ct-r)  . 

The  system  parameters  identified  Include  Cf  ,  C  and  g.  aa  above, 
and  the  additional  parameter  g2»  The  time-delay,  r,  was  fixed  at 
r  -  .05  sec.  This  selection  was  made  based  on  a  very  preliminary  scan 
of  the  data.  In  this  case  the  requrled  initial  data  Is  the  value  . 
Xj(0)  end  the  Initial  history  of  x^(*)  on  the  Interval  [-r,0] .  gather 
than  allow  the  computer  algorithm  complete  freedom  of  choice  in 
determining  the  initial  data  several  'experiments'  were  performed  with 
operator  selected  values.  The  values  selected  ware  x^(0)  -36., 
XjC-r/N)  -  4.0,  Xj(-2r/H)  -  .25  and  all  other  Xj  'knot'  values  are 
sero.  In  the  results  reported  the  SPLINE  REDUCED  procedure  wss  used 


TABLE  6.3.  Parameter  estimates  for  the  hereditary  model;  Example  6.2. 


iteration 

c>. 

C 

*1 

*2 

0 

6.28 

.500 

-12.2 

0 

1 

11.72 

.739 

-70.0 

13.7 

5 

9.36 

1.512 

-270.5 

-113.5 

10 

9.07 

1.658 

-309.8 

-104.6 

20 

8.81 

1.496 

-274.8 

-62.5 

49 

8.77 

1.412 

-255.1 

-47.5 

90 

8.76 

1.411 

-254.8 

-47.8 

Results  from  the  parameter  identification  algorithm  for  this 
case  are  shown  In  Table  6.3  and  in  Figure  6.3.  The  startup  fits 
are  the  same  as  In  the  previous  example  and  are  not  repeated.  While 
the  data  fits  for  the  hereditary  model  exhibit  some  Improvement  over 
the  previous  example  the  quality  of  the  matches  is  still  rather  poor. 
As  noted  above  It  seems  that  one  or  more  of  the  'fixed'  parameters 
has  been  Incorrectly  specified. 
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CONCLUSIONS 

The  state  space  model  (43)  offers  a  generalization  of  earlier 
formulations  of  dynamic  models  for  the  system  shown  In  Figure  1.1 
[13,  16,  17].  The  Inclusion  of  a  hereditary  dynamic  model  for  the 
'aerodynamic'  states  Is  significant  in  that  several  recent  Investigations 
have  shown  that  fundamental  unsteady  aerodynamic  analysis  produces 
precisely  such  hereditary  models  [1,  2,  9]. 

In  order  to  deal  efficiently  with  the  problem  of  parameter  estimation 
for  the  model  (43)  some  earlier  numerical  algorithms  for  identification 
of  hereditary  systems  [4,  5,  8]  were  reworked  to  Incorporate  separation 
of  'delayed'  and  'non-delayed*  states.  This  refinement  results  in 
significant  savings  In  computer  costs  and  increases  the  scope  of  pro¬ 
blems  for  which  the  methods  are  easily  Implemented. 

Some  results  from  preliminary  studies  of  identification  using 
actual  wind-tunnel  data  indicate  some  improvement  when  hereditary 
modelling  is  used.  These  results  must  be  viewed  with  caution  until 
more  systematic  studies  are  performed. 
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